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Preface 


This  Final  Report  contains  the  M.S.  Thesis  of  Mr.  Douglas  Hlmberger, 

/ 

submitted  for  publication  December  1979.  -4The  work  describes  the  mathematical 
background  and  the  physical  principles  which  were  used  to  perform  calculations 
of  sonic  reflectivity  from  the  water  side  of  a  water-ice-air  layered  system 
where  the  angle  of  incidence  is  varied  to  produce  changes  of  the  magnitude 
and  Che  phase  of  the  incident  sonic  beam  as  it  is  reflected  from  the  ice. 

It  is  shown  that  evaluation  of  these  changes  in  reflectivity  can  be  used  to 
determine  the  thickness  of  the  ice  without  having  to  be  in  direct  mechanical 
contact  with  the  ice. 

The  computation  of  modulus  and  phase  of  the  reflected  beam  is  accomplished 
via  a  modified  Simplex  Computation  Method.  The  basic  concepts  of  the  method 
are  describe^  as  are  the  modifications  incorporated  to  make  the  program 
useful  for  calculations  of  ice  thicknesses  which  one  would  ordinarily  expect 
to  encounter  in  the  Arctic  Ocean,  -r  *v,  *  •  /*;  >» 

The  Report  shows  that  there  are  ranges  of  angle  of  sonic  Incidence  and 
ranges  of  the  expected  product  ice  thickness  times  sonic  frequency  where  the 
results  of  reflectivity  measurements  are  sufficiently  characteristic  so  that 
variations  in  the  mechanical  properties  of  the  ice  (as  for  instance,  longi¬ 
tudinal  and  shear  wave  propagation  velocities)  influence  the  interpretation 
of  the  results  to  a  negligible  extent  so  that  an  assessment  of  the  ice  thick¬ 
ness  is  in  principle  possible  for  any  kind  of  ice.  It  is  also  pointed  out 
In  which  ranges  of  angle  of  incidence  and  sonic  frequency  times  expected 
ice  thickness  the  results  of  reflectivity  measurements  will  depend  highly 
on  a  previous  knowledge  of  the  mechanical  properties  of  the  ice,  i.e.,  which 
ranges  should  be  avoided. 


Examples  of  calculations  and  some  experimental  measurement  points  are 
given  to  illustrate  the  usefulness  of  the  technique.  A  complete  computer 


program  is  given  which  can  be  used  to  calculate  any  combination  of  parameters 
of  interest.  A  number  of  figures  are  presented  to  illustrate  the  influence  of 
changes  in  the  parameters  on  the  relfectivlty  results  and  thus  on  the  degree 
of  reliability  of  thickness  measurements. 

W  - 

Walter  G.  Mayer 
Principal  Investigator 


Washington,  D.C.,  December  1979 
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CHAPTER  1 


n 


INTRODUCTION 

Knowledge  of  the  propagation  characteristics  of  sound  waves, 
particularly  of  those  in  the  ultrasonic  region,  is  critical  in  many 
areas  of  scientific  research.  Notable  fields  among  these  are  seismic 
studies,  nondestructive  testing,  and  sonar  detection.  An  evaluation 
of  sound  reflection  from  and/or  trananission  through  a  particular  me¬ 
dium  may  reveal  information  about  a  number  of  physical  parameters  of 
that  medium.  One  of  these  parameters  is  the  thickness  of  a  solid  re¬ 
flector  bounded  by  a  liquid.  Firestone  [1]  observed  that  a  sound  wave 
impinging  on  a  flat  solid  plate,  inner sed  in  a  liquid,  will  be  reflected 
or  transnitted  depending  on  the  angle  of  incidence.  Maximum  trans¬ 
mission  will  occur  if  the  incident  sound  wave  excites  a  normal  mode 
of  vibration  of  the  plate.  These  normal  modes  of  vibration  were  de¬ 
scribed  by  Lamb  [2]  who  found  that  the  number  of  possible  modes  and 
their  propagation  velocities  are  determined  by  the  elastic  properties 
of  the  solid,  the  frequency  of  excitation,  f,  and  the  thickness  of  the 
plate,  d.  Experimentally  [3] ,  the  mode  structure  does  not  change  sig¬ 
nificantly  if  cne  considers  a  solid  plate  vibrating  in  a  vacuum  (the 
approach  used  by  Lamb)  or  in  a  fluid  as  long  as  the  density  of  the 
solid  is  at  least  three  or  four  times  greater  them  the  density  of  the 
surrounding  medium. 

j_  Huang  [4]  developed  a  set  of  formulas  which  shew  that  the  mode  __ 


IN^wiw  u.***w». 


-  2  - 


^structure  changes  significantly  when  the  solid  plate  is  bounded  by  a 
liquid  on  one  side  and  air  on  the  other.  Of  particular  interest  is  the 
case  where  such  an  asyimetric  loading  exists  and  where  the  density  of 
the  solid  is  less  than  the  density  of  the  liquid.  An  example  of  such 
a  system  is  an  ice  sheet  floating  on  water.  Since  the  ice  sheet  exhi¬ 
bits  normal  inodes  of  vibration  when  excited  by  an  underwater  acoustic 
wave,  it  is  conceivable  that  one  can  determine  the  thickness  of  the  ice 
from  the  results  of  reflectivity  measurements  at  various  angles  of  in¬ 
cidence. 

This  thesis  presents  a  method  of  determining  computationally  the 
reflectivity  as  a  function  of  angle  of  incidence  for  a  three  layer, 
three  substance  system,  specifically  for  the  air/ice/water  subsystem. 
Calculations  are  presented  for  an  fd  (frequency-thickness  product) 
range  up  to  4  x  106  Hz*nm. 
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CHAPTER  II 


1 


BACKGROUND 

The  basic  starting  point  for  the  study  of  sound  reflectivity  char¬ 
acteristics  is  the  analysis  of  the  behavior  of  sound  impinging  on  a 
flat  interface  between  two  differing  substances.  In  general,  a  sound 
wave  (or  beam)  incident  at  a  flat  interface  obeys  the  law  of  specular 
reflection;  that  is,  the  angle  of  reflection  is  equal  to  the  angle  of 
incidence.  However,  there  are  several  instances  [5]  where  this  does  not 
hold  true.  In  addition  to  the  wave  being  reflected  specularly  (and 
being  transmitted  to  sane  degree) ,  a  surface  wave  may  be  set  up  in  the 
solid  far  several  angles  of  incidence  (assuming  one  of  the  layers  is  a 
solid) .  These  are  waves  that  propagate  along  the  interface  between  the 
solid  and  the  liquid.  These  "leaky"  surface  wave.-'  are  attenuated  as 
they  re-radiate  energy  back  into  the  liquid.  The  generation  of  leaky 
surface  waves  at  these  special  angles  of  incidence  gives  rise  to  an  un¬ 
usual  type  of  reflection,  known  as  "non-specular " ,  and  gives  an  indi¬ 
cation  that  the  other  parameters  of  the  system  (thickness,  frequency) 
are  in  a  particular  arrangement.  To  determine  when  these  special  re¬ 
flections  occur,  one  must  first  find  a  functional  description  of  the 
reflection  and  solve  for  zero  reflection  (maximum  transmission) . 

The  most  simple  system,  which  Lord  Rayleigh  [6]  studied,  is  one  of 


an  infinitely  thick  solid  plate  bounded  by  a  vacuum.  He  determined 
that  the  solid  can  support  one  surface  wave,  the  Rayleigh  vave.  Hew- 
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Ever,  this  surface  wave  cannot  be  excited  by  an  incident  bean  as  there  ~ ^ 
is  no  median  to  propagate  such  a  sound  beam  to  the  solid. 

Because  it  was  shown  to  have  similar  vibrational  modes  [3] ,  the 
next  system  to  consider  is  that  of  an  infinitely  thick  plate  bounded  by 
sane  liquid.  Schoch  [5]  and  Brekhovskikh  [7]  studied  sonic  reflection 
at  the  interface  foamed  by  this  liquid/solid  systan  (or  I/S  system)  and 
reported  non- specular  reflections.  Figure  1  shows  the  law  of  specular 
reflection  far  the  incident  beam  and  the  transmission  of  part  of  the 
beam  by  both  shear  and  longitudinal  waves  in  the  solid.  The  directions 
of  propagation  of  the  various  waves  can  be  found  fran  Snell's  law, 
given  by 


sxn  0 


sin  e. 


sin  0_ 


where  i  denotes  liquid,  s  denotes  shear,  and  d  denotes  longitudinal. 

For  most  liquid/solid  systems  the  magnitudes  of  the  velocities  are  re¬ 
lated  by  the  inequality  v^  >  vg  >  v^  .  Using  these  facts  and  the  for¬ 
malism  of  Pitts  [8] ,  one  can  arrive  at  a  reflection  coefficient  (  a 
ratio  of  the  amplitudes  of  the  reflected  wave  and  the  incident  wave)  as 
a  function  of  the  angle  of  incidence  and  the  densities  of  the  liquid 
and  the  solid.  The  reflection  coefficient,  as  formulated  by  Pitts,  is 
given  fcy 
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Figure  1  —  Reflection  and  transmission  of  sound  wave 
at  a  liquid/solid  interface* 
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R(kx) 


=  (ll3  ' 

<* l  -  ax»2  +  *  ‘■“sV" 


where  the  following  definitions  are  used 


k  =  wave  number, 


K  =  ktcos  0  =  k  vL-sin  0  , 


K  =  k  cos  0  =  k  ' (v./v_) A  -  sin  9  ,  1 

3  S  •  S 

Kj=  k^cos  9  =  k  /(vyvj)^  -  sin^  0  ,  ^ 

0  =  angle  of  incidence  , 

p  =  density  of  liquid/  density  of  solid.  ^ 

Equation  (2)  has  a  zero  for  a  specific  angle  of  incidence  (in  a 
system  where  all  other  parameters  are  defined) .  It  is  at  this  angle 
where  a  surface  wave,  the  leaky  Rayleigh,  can  be  excited.  There  are, 
however,  certain  systems,  as  found  by  Brewer,  Himberger,  and  Mayer  [9] , 
where  this  solution  may  not  exist  due  to  certain  combinations  of  v(  , 
vg  ,  Vj  ,  and  p.  For  the  systems  where  a  solution  does  exist,  a  simple 
scan  or  search  method  would  be  sufficient  to  solve  for  the  zero.  The 
function  is  also  simple  enough  that  a  type  of  gradient  or  derivative 
method  could  be  used  (this  type  of  method's  use  depends  on  the  ease 
•with  which  the  function '  s  derivative  can  be  found) . 

The  systsn  becomes  more  complicated  when  the  solid  is  made  to  be 
|of  finite  thickness  and  is  loaded  on  both  sides  by  a  liquid,  in  which  _j 


r  ~i 

case  two  boundaries  must  now  be  considered.  Figure  2,  taken  fran  Pitts 
[8]  ,  shows  this  situation  in  terns  of  the  previously  defined  quanti¬ 
ties.  The  coefficient  formulated  by  Pitts  [8]  for  this  system  is  given 
by 

R(k  )  =  N/(f  f  )  (7) 

A  S  a 


where 


f 


s 


(k2  -  2k2) 2  (1  +  °°SP)  +  4k2K_K 

'  S  X  -  X  s  < 

sin  P 


(1  +  cosQ) 
sin  Q 


-  i'Vd 


K 


(8) 


f.  =  <k?-2k2)2  d-a&.K,  a-°°sQ) 

3.  s  x  —  x  s  a  . 

sin  P  sin  Q 

-  «kiKd 
___________  9 

K 

N  =  (kj  -  2k2)4  +  16k£^K2  .  0Vk^K2 

+  8(k2  -  2k2)  Vn,  (1-cosPcosQ)  /  (sinPsinQ) , 

S  X  X  S  Cl 

P  =  d  kj  cos  , 

Q  =  d  k  cos  ©  ,  and 

s  s 


(9) 


(10) 

(11) 

(12) 


d  =  thickness  of  plate. 
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The  relationship  between  the  angle  of  incidence,  Q,  and  the  thickness 
of  the  solid  plate  is  canplicated.  This  situation  requires  a  more  so¬ 
phisticated  type  of  optimizing  technique. 

One  finds  several  zeroes  for  (7)  for  a  given  fd  or  angle  of  inci¬ 
dence.  These  solutions  correspond  to  the  modes  of  vibration  the  inci¬ 
dent  beam  is  capable  of  exciting  in  the  solid  plate.  An  example  system 
is  the  water/brass /water  system  of  Pitts  [8] .  Figures  3  and  4  show  the 

real  part  of  the  zeroes  of  the  f  and  f  terms  for  different  fd's  and 

s  a 

incident  angles.  Although  solutions  to  these  terms,  equations  (8)  and 
(9) ,  represent  solutions  to  the  denaninator  of  the  reflection  coeffi¬ 
cient,  and  thus  are  poles,  it  has  been  shewn  [8]  that  the  solutions  to 
the  numerator  (the  zeroes)  have  the  same  real  sin  0  value.  It  has  also 
been  shown  [8]  that  these  pole-zero  pairs  are  complex  conjugates  of 
each  other  when  the  imaginary  part  is  also  included.  It  is,  however, 
the  real  part  of  the  angle  that  will  be  measured  when  a  non-specular 
reflection  is  observed.  The  curves  in  figs.  3  and  4,  taken  from  Pitts 
[8]  ,  are  referred  to  as  modes,  representing  the  modes  of  vibration  of 
the  solid  plate. 

The  relation  given  by  (7)  becomes  more  canplicated  when  one  goes 
to  an  even  more  general  system  in  which  the  plate  is  surrounded  by  two 
dissimilar  fluids,  or  by  a  fluid  and  a  gas,  as  in  the  air/ice/water 
case.  There  is  no  distinction  between  "syrnnetric"  and  "anti-synmetric" 
solutions  to  the  reflection  coefficient  (Pitts'  solutions  to  the  fg 
gnd  f  terms  respectively)  because  neither  the  system  nor  the  coeffi-  ^ 
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'cient  are  symmetric.  Cne  can  no  longer  hreak  up  the  coefficient  into 
two  siirple  terms  and  find  the  zeroes  for  each  of  these  terms.  The  en¬ 
tire  equation  nust  be  solved  simultaneously.  Fran  Huang  [4] ,  the  co¬ 
efficient  for  this  type  of  system  is  given  by 


R  =  N/D  , 

where 

N  =  A(B-C)  +  (E+F)G  , 

D  =  A(B+C)  -  (E-F)G  , 

and  -  0  . 

A  =  2  {  a  -  [F(sin  {  irfd/v^  [  (v^/v^)  ^  -  sin  eP  })  • 


(13) 

(14) 

(15) 


(sin(  tt  fd/v1  [  (v1/vs) 2  -  sin2  ©]  ^  })]}  ,  (16) 

B  ={  1-2 [  (vg/v1) sine]2}2  sin(  wfdA^ t  (v1/v£j) 2  -  sin2e]'5}  • 
cos  {  irfd/v^  [  (v^) 2  -  sin2e]  * }  + 

{  4 [  (Vg/v^ 4sin2e]  [  (v^/Vj)  2-sin20]  ** [  (v-j/Vg)  2-sin2e] ^ }  • 

cos  { Trfd/v^ [  (Vj/Vj)  2-sin20] ** }  sin  {  Trfd/v^ [  (v^/Vg)  2-sin2e] , 

(17) 

C  =  E  cos  {irfd/v1[(v1/vd)2-sin2e]Js}cos  { irfdA^  [  (v-^/Vg)  2-sin20]  S  , 

(18) 

E  =  i  (P1/P2)  (  [  (Vj/Vj)  2-sin20]  V (1-sin2 ©)**}  ,  (19) 

F  =  i  (p3/p2)  {  [(v^J^s^elVUv^J^sin2^15}  ,  (20) 


L 
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G  =  {1-2 [  C^g/v^ sin  9} 2 }2  sin{ irfc3/v1  [  (Vj/vg)  2-sin20]  **}  • 
cos  {irfd/v1[(v1/vs)2-sin20]is}  + 

4 [  (vg/v1)  4sin20]  [  (Vj/Vj)  2-sin20] ** [  (v^Vg) 2-sin20] ^  • 
sin  { irf d/v^  [  ( v^/v^) 2— sin20  3  ^ }  cos  { Ttfd/v^  [  (v-^/v^)  2-sin20] 

(21) 

n  =  {l-2[(vs/v1)sin  0]2}2  cos  { irfd/v.^ [ (v^/v^) 2-sin20] ^ }  • 

sin{  Trfd/v1  [  (v1/vs)  2-sin20]  ** }  + 

4 [  (vg/v]L)  4sin20]  [  (^/Vj)  2-sin20] ** [ (v^/Vg)  2-sin29] **  ,  (22) 

=  velocity  in  medium  1  (liquid  or  gas)  , 

-  velocity  in  medium  3  (liquid  or  gas)  , 
pi  =  density  of  medium  1  , 

p2  =  density  of  medium  2  (solid)  , 

P3  =  density  of  medium  3  . 

Equations  (16)  through  (22)  are  expressed  on  terms  of  velocities 
rather  than  wave  numbers  as  in  (2),  (8)- (12) .  Optimizing  (13)  re¬ 
quires  more  than  a  simple  technique. 

The  air/ice/water  system  is  a  special  case  of  an  L/S/L  system 

because  v  and  v  may  be  related  by 
1  s 

L 
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vg  >  ,  (in  most  cases) 

vg  <  v  ,  (in  rare  cases) 

where  again  v^  refers  to  the  sound  velocity  in  the  water  and  vg  refers 
to  the  shear  velocity  in  the  ice.  Clearly,  the  sound  velocity  in  air 
is  always  less  than  v  .  Moreover,  this  system  differs  fran  most  VS/L 
sy stons  because  the  density  of  the  solid  is  less  than  that  of  the  load¬ 
ing  liquid  (pice  =0.917,  =1.02). 

Although  the  shear  wave  velocity  in  the  ice  may  be  less  than  the 
sound  velocity  in  water  [10] ,  the  present  paper  is  concerned  only  with 
L/S/L  sy stans  where  vg  ^  v^  ,  the  more  common  situation. 

A  computational  method,  the  simplex,  was  adapted  to  solve  (13) 
representing  this  complicated  general  system.  The  interpretation  of 
the  results  then  determined  the  feasibility  of  finding  the  thickness 
of  an  ice  layer  by  observing  the  variations  in  the  reflection  as  a 
function  of  sonic  frequency  and  angle  of  incidence. 
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CHAPTER  III 


THE  SIMPLEX 


To  solve  the  reflection  coefficient  expressed  by  (13)  in  the  most 
efficient  and  flacible  manner,  five  methods  were  investigated:  the 
Newton-Raphson  gradient  method  [11] ,  the  Rosenbrock  direct-search  me¬ 
thod  [12],  the  Powell  'sum  of  squared  residuals'  method  [13],  the 
Powell  direct  search  method  [14] ,  and  the  Nelder  and  Mead  simplex  dir¬ 
ect  search  method  [15] . 

One  of  these  minimization  techniques  was  excluded  from  consider¬ 
ation  as  a  result  of  initial  investigation .  The  Newton-Raphson  (and 
others  similar  to  it,  such  as  the  Fletcher-Pcwell  method  [16] )  is  a 
gradient  method  which  requires  that  not  only  the  function  itself  be 
defined  but  also  the  derivative  of  the  function  be  known.  For  many 
functions  this  is  not  a  problem,  but  upon  close  inspection  one  sees 
that  (13)  does  not  appear  to  lend  itself  to  this  type  of  technique  be¬ 
cause  of  the  complexity  of  the  function  and  the  difficulty  of  finding 
its  derivative. 

Powell's  'squared  residual'  optimization  technique  (and  others 
with  similar  characteristics  such  as  Barnes'  method  [17])  does  not  re¬ 
quire  that  the  derivative  of  the  function  be  known.  While  this  is  an 
advantage  in  many  cases,  Powell's  method  is  not  suited  for  all  func¬ 
tions.  Bax  [18]  states,  "It  may  not  be  possible  to  reformulate  every 
jCptimization  problem  as  the  solution  of  a  set  of  simultaneous  equa-  j 


-  15  - 


.1 


1 


tions  [which  Fcwell ' s  method  requires  for  the  ' squared  residuals ' ] .  .  . 
Certainly  it  is  often  not  feasible  to  attempt  this."  It  is  unlikely 
that  (13)  oould  easily  be  separated  into  the  desired  set  of  equations. 
For  a  complicated  function,  Bax  [18]  states,  "It  has  been  known  for 
the  confutation  of  all  residuals  for  a  single  set  of  parameter  values 
to  take  if)  to  1,000  times  as  long  as  the  organization  of  the  search." 
Even  if  the  set  of  equations  necessary  for  Powell's  method  was  avail¬ 
able,  the  factor  of  computation  time  would  ranove  this  type  of  tech¬ 
nique  fran  consideration. 

Sane  type  of  'direct  search'  method  is  more  likely  to  minimize  (13) 
without  the  problems  associated  with  the  previous  methods  discussed. 

The  direct  search  methods  investigated  were  Rosenhrock's  [12],  Powell's 
[14] ,  and  Nelder  and  Mead's  [15] . 

Kcwalik  and  Osbcume  [19]  compare  Nelder  and  Mead's  method,  the 
simplex,  to  Rosenbrock's  method  and  state  the  two  methods  have  "can- 
parable  efficiency  when  tested  on  problems  with  a  small  number  of  in¬ 
dependent  variables.  .  .  however,  .  .  . [for  sane  specific  functions] 
the  Simplex  method  has  shown  superiority  over  [Rosenbrock's  method]". 
Kcwalik  and  Osbourne  also  agree  with  Nelder  and  Mead  in  the  latters ' 
evaluation  that  the  simplex  method  held  an  advantage  over  Powell's 
method  in  having  a  faster  convergence  for  several  example  functions. 

In  addition,  should  the  function  have  two  or  more  solutions  within 
close  proximity,  Nelder  and  Mead  state  that  the  simplex  "will  con- 
jverge  even  when  the  initial  simplex  straddles  two  or  more  valleys, 
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a  property  which  is  not  shared  by,  e.g. ,  Powell's  method."  Because 
the  forms  of  the  modes  (the  solutions)  of  (13)  are  not  known,  the  abil¬ 
ity  of  the  method  to  converge  regardless  of  the  initial  guess  is  use¬ 
ful.  The  simplex  method,  in  Nelder  and  Mead's  words,  is  "highly 
opportunist,  in  that  the  least  possible  information  is  used  at  each 
stage  and  no  account  is  kept  of  past  positions.  No  assumptions  are 
made  about  the  surface  except  that  it  is  continuous  and  has  a  unique 
minimum  in  the  area  of  the  search."  In  addition,  the  simplex  method, 
according  to  Nelder  and  Mead,  is  " computationally  compact".  For  all 
of  the  above  reasons,  the  simplex  was  chosen  as  the  minimizing  tech¬ 
nique  to  solve  for  the  zeroes  and  poles  of  (13) . 

The  simplex  method  of  optimization  was  developed  by  Himsworth, 
Spendley  and  Hext  [20  ] ;  it  was  later  studied  in  detail  by  Nelder  and 
Mead  [15] .  Although  the  technique  was  originally  developed  for  use 
in  business  and  economics  (plant  management,  etc.) ,  the  method  was 
used  in  mathematics  and  science  fields  to  seme  degree  in  the  early 
1960's.  Its  use  in  technical  fields  has  decreased  however,  and  to¬ 
day  the  method  is  utilized  almost  exclusively  by  the  business  world. 

The  sinplex  is  a  'steep  ascent'  method  that  needs  no  derivatives 
and  simply  forms  a  geometric  figure,  the  simplex  (from  which  its  name 
is  derived) ,  and  calculates  the  actual  functional  value  at  the  ver¬ 
tices  of  this  figure.  The  definition  of  a  sinplex  is  given  below  as 
a  direct  quote  from  Kcwalik  and  Osbourne  D-9  ]  ? 
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.  .a  set  of  n+1  points  in  n-dimensional  space  forms  a  simplex.  When 


these  points  are  equidistant  the  simplex  is  said  to  be  regular."  Ihe 
role  of  the  simplex  is  to  find  the  minimum  of  a  particular  function  by 
"sliding"  along  the  functional  value,  whether  this  sliding  is  done 
along  a  one-dimensional  line,  a  two  dimensional  plane,  or  a  three  di¬ 
mensional  surface.  The  simplex,  as  stated  by  Nelder  and  Mead  tL5  ] , 
"adepts  itself  to  the  local  landscape,  elongating  down  long  inclined 
planes,  changing  direction  on  encountering  a  valley  at  an  angle,  and 
contracting  in  the  neighborhood  of  a  nunimum. "  For  added  flexibility 
in  finding  minima,  Nelder  and  Mead  generalized  the  simplex  to  make  it 
non-regular,  that  is,  not  necessarily  regular  or  symmetric  in  nature. 

The  algorithmic  process  of  the  simplex  is  given  in  fig.  5,  taken 
from  Nelder  and  Mead  D-5  ] .  The  simplex  has  three  basic  operations; 
reflection,  expansion,  and  contraction.  The  definitions  of  these  op¬ 
erations  and  of  basic  terms  are  given  by  Kcwalik  and  Osbourne  0.9  ]  as 
follows; 

"(1)  is  the  vertex  which  corresponds  to  f (x^)= 
max  f(x^),  where  i=l,2, . . . ,n+l. 

(2)  x  is  the  vertex  which  corresponds  to  f(x  )= 

s  s 

max  f (x^) ,  where  i  /  h. 

(3)  x„  is  the  vertex  corresponding  to  f (x^)= 
min  f(x^),  where  i=l,2,...,n+l. 

(4)  xQ  is  the  centroid  of  all  x^,  i  ^  h  and  is 
given  by 
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We  now  define  the  three  basic  operations  used 
in  the  method: 

(1)  Reflection,  where  x^  is  replaced  by 

xr  =  (1  +  a)XQ  -  a  x^ 

where  the  reflection  coefficient  a  >  0  is  equal 

to  the  ratio  of  the  distance  [x  x  ]  to  [x.  x  ]  . 

r  o  no 

(2)  Expansion,  where  xr  is  expanded  in  the  di¬ 
rection  along  which  a  further  improvanent  of  the 
function  value  is  expected.  We  use  the  relation 

xe  =  yxr  +  (1  -  y)xq, 

where  the  expansion  coefficient  y>  1  is  the  ratio 
of  the  distance  lxexQ]  to  [xrxQ]  • 

(3)  Contraction,  by  which  we  contract  the  simplex, 

xc  -  0  \  +  (I"  S>x0' 

where  the  contraction  coefficient  S  is  the  ratio  of 

the  distance  [x  x  1  to  [x,  x  1  and  satisfies 
co  no 

0  <  B  <  1. 

As  we  have  mentioned,  the  method  can  be  viewed  as 
the  moving,  shrinking,  and  expanding  progress  of  the 
simplex  toward  the  minimum.  This  motion  is  accam- 
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plished  in  the  following  way: 

(i)  An  initial  simplex  is  formed,  and  the  func¬ 
tion  is  evaluated  at  each  of  the  vertices  in  order 

to  determine  x,  ,x  ,x  ,  and  x  . 

n  s  i  o 

(ii) We  first  try  reflection  and  evaluate  the 
function  at  the  reflected  point. 

(iii)  If  f(xg)  >_  f(xr)  >_  f  (x^) ,  then  we  replace 
x^  by  xr  and  restart  the  process  with  the  newly 
formed  simplex. 

(iv)  However,  if  f(xr)  <  f(x^),  we  may  expect  that 

the  direction  x  -x  could  give  us  an  even  lower 
r  o 

value  of  the  function  if  we  move  further.  There¬ 
fore  we  expand  our  new  simplex  in  this  direction. 

The  expansion  succeeds  is  f(xj  >  f(x  ),  and  in 

36  6 

this  case  x^  is  replaced  by  xr,  and  in  either  case 
we  restart  the  process  fran  our  new  simplex. 

(v)  if  the  reflection  move  (ii)  yields  xr  such  that 
f  (x^)  >  f  (xr)  >  f  (xg) ,  we  replace  x^  by  xr  and 
make  the  contracting  move.  This  replacement  is 

not  executed  when  f (xr)  >  f (x^) .  After  the  contrac¬ 
ting  move  we  compare  the  f (x^)  and  f (xc) .  If  f (x^) > 
f (xc) ,  we  consider  that  the  contraction  is  success¬ 
ful,  x^  is  replaced  by  xc,  and  we  start  from  the 
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new  simplex. 

In  a  case  of  failure,  i,e.,  f(x^)  <  f(xc),  the  last 
simplex  is  shrunk  about  the  point  of  the  lowest  func¬ 
tion  value  x  by  the  relation 

xi =  ls(xi +  xt  > 

and  we  begin  from  (i) . 

The  stopping  criterion  suggested  by  Nelder  and 
Mead  is 

i=l 

II 

where  e  is  seme  anall  preset  number. 


n+1 


{  1/n 


[f(xi)  -  f(xo)]‘ 
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CHAPTER  IV 


MODIFICATIONS  CF  THE  SIMPLEX 

The  simplex  routine  was  executed  in  a  FORTRAN  IV  computer  pro¬ 
gram  for  use  on  both  the  IBM  370/148  canputer  at  the  Academic  Compu¬ 
tation  Center  at  Georgetown  University  and  a  DEC  PDP-MINC  11  system 
at  Harry  Diamond  Laboratories.  FORTRAN  was  used  because  of  its 
capacity  for  handling  involved  calculations  with  ease  as  well  as 
its  ability  to  deal  with  complex  numbers.  This  was  necessary  due  to 
the  fact  that  (13)  has  solutions  that  are  in  general  complex. 

The  program  itself  was  based  on  a  simplex  program  developed  by 
Huang  [4]  to  deal  with  similar  problems;  both  are  derived  from  Nel- 
der  and  Mead's  flowchart  as  given  in  chapter  III  of  this  paper. 

The  functional  changes  made  in  Huang's  program  involve  the  fol¬ 
lowing  parts  of  the  routine; 

(i)  When  no  'new  iminimum'  is  reached  by  a  reflection  (the 
operations  identified  in  section  (v)  of  Kcwalik  and 
Osbourne,  as  quoted  in  chapter  III) ,  the  new  contrac¬ 
ted  simplex  parameters  are  specified. 

(ii)  When  setting  up  the  initial  simplex,  the  centroid 
calculation  is  set  as  specified  by  Nelder  and  Mead 
[153. 


These  changes  in  Huang '  s  program  made  the  new  simplex  routine  both 
jjxrnputationally  faster  and  numerically  more  consistent.  The  number 


of  iterations  needed  to  reach  the  minimum  of  a  test  function  was  re¬ 


duced  fran  an  average  of  150  to  80.  The  minimum  found  in  the  test 
function  was  consistent  to  five  decimal  places  even  when  the  iter¬ 
ation  was  begun  at  several  different  starting  points. 

Nelder  and  Mead  [15]  state  that  a  "general  problem  encountered  by 
all  minimization  methods  is  that  of  false  convergence  at  a  point  other 
than  the  minimum.  This  difficulty  has  been  found  using  the  simplex 
method...".  It  is  because  of  these  'false  minima'  that  the  ability  to 
trace  a  mode,  such  as  one  similar  to  those  in  figs.  3  and  4,  becomes 
important.  Lacking  this  ability,  it  would  be  difficult  to  differen¬ 
tiate  between  local  minima  and  the  physically  meaningful  solutions. 
Several  features  were  added  to  the  simplex  which  remove  this  difficul¬ 
ty.  The  first  was  a  set  of  signal  'flags'  telling  the  investigator 
what  series  of  reflections,  attractions ,  etc.  of  the  simplex  were 
used  to  reach  a  final  minimum.  This  was  done  in  order  to  predict  to 
which  minimum  the  simplex  would  iterate  from  a  given  starting  point. 

Secondly,  the  program  was  put  into  two  separate  forms;  a  concise, 
simple  'scan'  of  a  given  fd  for  all  values  of  sin  0,  and  a  rigorous 
mode-tracing  form  that  not  only  followed  a  particular  mode  but  gave 
information  about  adjacent  modes.  The  first  form  was  used  to  find  the 
starting  points  of  any  and  all  modes  for  a  set  of  parameters  and  the 
second  was  used  to  generate  full  sets  of  curves  for  these  parameters. 

Several  special  features  were  added  to  the  mode-tracing  form. 


^Though  the  simplex  inherently  samples  the  function's  three  dimensional^. 
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space  (in  the  present  case  these  dimensions  being  the  product  fd, 
the  value  of  sin  0,  and  the  functional  value  for  any  of  these  two 
parameters)  until  it  finds  a  valley  into  which  to  go,  the  new  form  of 
the  simplex  was  'farced'  to  search  in  a  given  volume  of  a  suspected 
minimum  even  if  another  mode  was  nearby.  While  this  search  was  in 
progress,  other  minima  found  (or  even  tended  towards)  were  indicated 
in  the  output  so  that  a  better  knowledge  of  the  function  in  the  three 
dimensional  space  under  scrutiny  would  be  gained.  When  dealing  with 
the  beginning  of  the  mode  in  question,  this  range  of  forced  iteration 
was  determined  by  the  initial  fd  scan  form  of  the  simplex.  In  follow¬ 
ing  the  mode,  once  a  solution  was  found  in  the  sampling  volume,  the 
range  of  forced  iteration  was  determined  by  the  last  point  found  and  a 
judicious  choice  of  limits  set  by  the  operator  watching  the  tracing 
of  the  mode. 

When  the  simplex  could  not  iterate  to  a  minimum  within  the  given 
limit  of  iterations  (determined  by  the  amount  of  computer  time  avail¬ 
able)  ,  the  output  also  included  the  intermediate  value  reached  so  that 
sane  information  would  be  gained  about  a  troublesome  range  of  the  func¬ 
tion.  This  point  could  then  be  re-examined  after  slightly  varying  the 
starting  point,  the  iteration  step  size,  or  both. 

The  function  described  by  (13)  varies  in  such  a  way  as  to  make  the 
above  additions  to  the  simplex  useful.  For  most  test  functions  such 
as  those  used  by  Box  [18] ,  Kowalik  and  Osbourne  [19] ,  and  this  author 
|J during  the  preliminary  sinplex  testing) ,  no  surface  around  a  minimum  _j 


had  a  form  that  caused  the  simplex  method  problems  with  respect  to 
local  minima,  adjacent  modes,  etc.  For  the  function  given  by  (13)  how¬ 
ever,  this  was  not  true.  Figure  6  shows  two  different  example  minima, 
each  with  considerably  different  approach  slopes.  If  a  function  varied 
in  a  manner  shown  in  fig.  6(b) ,  a  small  change  in  the  independent  var¬ 
iable  (in  the  present  case,  the  fd  or  sin  ©)•  would  cause  a  great  change 
in  the  dependent  variable  (the  functional  value  in  the  present  case) . 
This  rapid  variance  would  cause  the  simplex  to  miss  or  'step  over'  one 
mode  in  favor  of  another  close  and  more  gradually  approached  mode,  such 
as  the  one  shown  in  fig.  6(a) .  When  this  happened  in  tracing  a  mode, 
a  physically  meaningful  solution  was  regarded  at  times  as  a  local  min¬ 
imum.  Unfortunately,  the  form  of  the  solution  to  the  reflection  coef¬ 
ficient  of  the  air/ice/water  case  has  many  such  steep  approaches  to 
minima.  An  example  of  this  'steep  valley'  situation  is  indicated  in 
Table  I  which  represents  a  zero  of  the  denominator  of  the  reflection 
coefficient  (eq.  (13))  when  the  system  has  a  density  ratio  of  the  first 
liquid  to  the  solid  of  0.3  (compared  to  1.1  for  normal  water/ ice  den¬ 
sity  ratio) ,  and  velocities  in  ranges  similar  to  the  normal  air/ice/ 
water  case. 


Table  I.  "Steep  Valley"  variation  of  function  value. 


Starting  sin  0 

0.540 

0.600 


Real  sin©  Imag.  sin©  Function  Val. 


0.562996 

0.562996 


-.0117405 

-.0117408 


-.568E-06 

-.293E-05 


v, =1500  v_=1550  v*=3500  v,=340  fd=0.65  p./p  =0.3 
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^lhe  approach  slope  is  so  steep  that  when  the  real  parts  of  sin  0  are 
equal  to  six  places,  a  change  in  the  sixth  decimal  place  of  the  imagi¬ 
nary  part  causes  a  change  in  the  functional  value  of  one  order  of  mag¬ 
nitude.  This  solution  was  missed  in  both  a  scan  and  a  mode  trace  when 
using  the  old  simplex  form. 

Besides  improving  the  reliability  of  the  simplex  in  terms  of  find¬ 
ing  solutions,  the  additions  enable  the  new  form  to  trace  complete  modes 
in  cne- tenth  the  time  of  the  old  form.  This  reduction  in  computation 
time  makes  the  simplex  a  reasonable  tool  for  solving  complicated  func¬ 
tions  such  as  (13) .  The  final  form  of  the  simplex  mode  trace  routine 
is  given  in  Appendix  A. 


L 


J 


CHAPTER  V 

RESULTS 

A.  TVP-substance  systan  results. 

The  modified  simplex  optimizing  technique,  as  described  in  chap¬ 
ter  4,  was  used  to  solve  (13)  for  parameters  of  a  water/brass/water 
systan,  previously  described  by  Pitts  [8] .  Table  II  contains  the  re¬ 
sults  of  the  simplex  and  those  of  Pitts. 


Table  II.  Water/Brass/Water  Reflection  Coefficient  Poles 

(fd  =  5) 


Pitts 

Simplex 

Re(sinG) 

im(sinO) 

Re (sine) 

Bn  (sine) 

0.002 

0.454 

0.290307E-02 

0.454711E00 

0.004 

0.302 

0.485337E-02 

0.302640E00 

0.182 

0.0075 

0.182645E00 

0.751577E-02 

0.273 

0.0017 

0.273108E00 

0.177726E-02 

0.342 

0.0001 

0.342847EOO 

0. 995920E-04 

0.408 

0.0048 

0. 408948E00 

0.487755E-02 

0.544 

0.0057 

0.544637E00 

0. 57863 9E-02 

0.656 

0.0034 

0. 656149E00 

0.347561E-02 

0.747 

0.0136 

0.747784E00 

0.136844E-01 

0.750 

0.0125 

0.750331E00 

0.125455E-01 

The  results  of  the  simplex  confirm  Pitts'  original  calculations  for 
the  water/brass/water  system.  The  results  also  confirm  that  (13)  re- 
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duces  to  (7)  when  the  same  liquid  bounds  the  plate  on  both  sides. 


B.  Three-substance  systgn  results. 

Values  for  air,  ice  and  water  parameters,  given  by  Sato  [10] ,  were 
substituted  in  (13)  and  this  function  was  run  in  the  simplex.  Figure 
7  shows  the  results  of  the  runs  for  the  air/ice/water  reflection  coef¬ 
ficient  for  fd's  of  zero  to  four.  The  graph  shows  a  general  farm  simi¬ 
lar  to  that  of  figs.  3  and  4  from  Pitts  [8] .  There  are  several  major 
differences,  however.  Two  of  the  modes  shown  in  fig.  7  are  not  con¬ 
tinuous  whereas  all  of  the  modes  in  previous  studies  have  been.  The 
modes  intercepting  the  abscissa  (representing  integer  half-wavelength 
"standing  wave"  solutions  for  incoming  beams  of  normal  incidence)  cor¬ 
respond  to  the  half-wavelengths  of  the  shear  wave  as  expected,  but  no 
intercepts  representing  the  longitudinal  wave  occur.  Finally,  one  of 
the  modes  has  a  definite  negative  slope,  suggesting  a  negative  group 
velocity.  Viktorov  [21]  anticipated  this  possibility,  but  it  has  not 
been  confirmed.  Due  to  the  above  differences,  several  changes  were 
investigated  in  an  effort  to  clarify  the  situation. 


C.  Mjustment  of  solid  density. 

Brower,  Himberger,  and  Mayer  [9]  showed  that  leaky  Rayleigh  sur¬ 
face  waves  cannot  be  generated  on  certain  systems,  specifically  the 
ice/water  system.  Zeroes  of  the  reflection  coefficient  are  known  to 
have  a  definite  relationship  to  ncn-specular  reflection  and  hence  to 
^surface  waves,  Since  the  ice/water  system  was  found  to  be  a  systen  j 
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Figure  7.  Mode  trace  for  air/ice/water  system.  Sato 
values  used}  v^=  3500  m/s,  vg=  1550  m/s, 

1500  m/s,  vair=  340  m/s,  fice=  .91 7, 
/water'  I'02'  /air'  °-001183- 


i  i'-i  A'  jir.fc 
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^where  the  leaky  Rayleigh  wave  cannot  exist,  a  similar  system  such  as 
the  air/ice/water  systen  could  well  have  differences  attributable  to 
this  unique  situation  which  does  not  exist  for  other  air/solid/liquid 
systems.  To  check  this  theory,  four  different  ice  densities  were  used 
in  (13),  giving  different  values  for  the  density  ratio  of  p  /p  .  Fig- 
ure  8  shows  the  "plateau"  from  Brewer  et.  al.  [9] ,  below  which  the 
leaky  Rayleigh  wave  exists,  above  which  it  does  not,  and  the  four  den¬ 
sity  ratios  —  two  above  and  two  belcw  the  "plateau".  This  change  of 
parameter  value  was  purely  of  investigative  nature  since  the  density- 
modified  "ice"  represented  not  ice  at  all,  but  a  solid  with  identical 
shear  and  longitudinal  wave  velocities  and  different  density. 

Several  changes  can  be  seen  in  figs.  9-12,  which  shew  the  solu¬ 
tions  to  (13)  for  the  four  density  ratios.  The  mode  that,  in  fig.  7, 
had  a  negative  slope  now  has  a  positive  slope  throughout  the  investi¬ 
gated  range.  This  same  mode  now  approaches  asymptotically  the  Ray¬ 
leigh  wave  velocity  for  this  systsn  as  predicted  by  Pitts  [8] .  The 
top  mode  which  was  discontinuous  in  fig.  7  is  new  part  of  this  new 
positive  slope  mode.  Sane  problems  are  still  present,  however.  One 
mode,  in  both  figs.  9  and  10,  has  a  negative  sloped  area  and  is  dis¬ 
continuous.  There  still  cure  no  modes  representing  the  longitudinal 
wave  intercepting  the  abscissa  (3in©=0) . 

D.  Adjustment  of  acoustic  impedance. 

The  acoustic  impedance  is  often  important  when  considering  re- 
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Figure  8.  "Plateau"  of  existence  for  Rayleigh  wave.  Points 

below  plateau  represent  systems  where  the  Rayleigh 
surface  wave  dpes  exist,  points  above  represent 
those  where  it  does  not.  Point  1  has  density  ratio 
of  0.1;  point  2,  0.3;  point  3,  1.11;  point  4,  11. 
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Mode  trace  for  air/ice/water  system  with  mod¬ 
ified,  ice  density.  Average  Sato  values  used, 
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is  below  "plateau"  (ratio=0.1)  where  leaky 
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Figure  10.  Mode  trace  for  air/ice/water  system  with  mod¬ 
ified  ice  density.  Average  Sat6  values  used, 
■Pice"  3-06,  Density  ratio  of  liquid  to  solid 

is  below  "plateau"  (ratio=0. 3)  where  leaky 
Rayleigh  exists. 
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FD  (FHEQIEMCT«TOIOOJE55) 

Figure  11.  Mode  trace  for  air/ice/water  system.  Average 
Sat6  values  used.  Density  ratio  of  liquid  to 
solid  is  above  "plateau"  (ratio=l.ll)  where 
leaky  Rayleigh  does  not  exist. 
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Figure  12.  Mode  trace  for  air/ice/water  system  with  mod¬ 
ified  ice  density.  Average  Sato  values  used, 

=  0.10.  Density  ratio  of  liquid  to  solid 


is  above  "plateau"  (ratio=ll)  where  leaky  j 

Rayleigh  does  not  exist.  I 


flections  from  an  interface.  A  matched-impedance  system,  or  one  in 

which  the  wave  velocity  times  the  density  of  each  substance  is  equal, 

\ 

generally  behaves  differently  than  one  that  is  not  matched.  Again,  a 
hypothesis  was  made  that  a  matched-impedance  systsn  between  the  ice  and 
water  would  eliminate  the  problem  of  di scon ti nous  and  negative  sloped 
inodes  as  shown  in  figs.  7  and  9-12.  To  consider  the  hypothetical  case 
of  matched  impecances  of  water  and  ice,  the  ice  density  was  again 
changed.  Figure  13  shows  the  solutions  to  the  impedance-matched  sys¬ 
tem,  the  normal  ice/water/air  system,  and  one  of  the  "belcw  plateau" 
density  ratio  systems.  The  impedance-matched  system  gives  results  that 
differ  little  frcm  the  normal  system.  Specifically,  the  previously 
discussed  problems  are  still  present,  leading  one  to  believe  that  the 
relatively  snail  difference  in  acoustic  impedances  between  ice  and 
water  is  not  the  reason  discontinuous  modes  exist. 


E.  Expansion  of  oontinuous-mode  range. 

The  curves  located  in  the  area  bounded  by  fd  of  1.75  and  4  and 
sinS  of  0.2  and  1.0  do  not  show  discontinuities  and  negative  slopes. 
Cue  to  the  frequencies  and  plate  thicknesses  ccmmonly  used  in  the 
laboratory,  this  area  of  solutions  is  the  area  most  often  verified 
by  experimentation.  For  these  reasons,  (13)  was  run  in  the  simplex 
using  high,  average,  and  low  shear  and  longitudinal  wave  velocities 
for  ice  (from  Satd  [10] )  in  this  bounded  area,  other  parameters 
(density,  etc.)  being  held  constant.  Shown  in  figs.  14  and  15  are  the 


Figure  13. 


Partial  mode  trace  for  air/ice/water  system 
using  three  different  ice  densities;  actual 
ice,  £*0.917  (solid  line);  impedance-matched 
ice,  j>=0. 987097  (dashed  line);  and  "below 
plateau"  ice,  f *3.06  (line  with  triangular 


1500  m/s,  vair=340  m/s.  Low  values  (dotted  lines)  are;  v^=34 00  m/s 
1500  m/s  (other  values  are  the  same). 
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Equality  of  the  imaginary  part  of  the  solution,  fig.  16  (a)  or  (b)  1 

may  each  be  considered  equally  correct  in  representing  the  physical  sit¬ 
uation. 

Pitts  [8]  has  also  stated  that  the  magnitude  of  the  imaginary  part 
of  the  solution  indicates  whether  the  zero  will  be  a  physical  as  well 
as  matherratical  solution  to  the  systan.  There  has  been  no  determina¬ 
tion,  however,  as  to  the  limits  of  the  magnitude  of  the  imaginary  part 
where  this  is  a  factor. 
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CHAPTER  VI 
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CCNCIUSICNS  AND  SUGGESTIONS  FOR 
FUTURE  RESEARCH 

Equation  (13)  represents  a  complicated  function  which  causes  prob¬ 
lems  with  optimization  not  found  in  simpler  functions  such  as  that 
represented  by  (7) .  Highly  varying  functioned  values,  portions  of 
inodes  with  negative  slope,  and  discontinuous  modes  all  contribute  to 
the  difficulties  experienced  in  solving  the  equations  for  air/ice/ 
water  system  inodes. 

The  simplex,  in  its  modified  form,  is  able  to  efficiently  optimize 
complicated  functions  such  as  (13).  With  operator  intervention  (watch¬ 
ing  imaginary  part  progression,  etc.) »  the  routine  can  trace  modes  as 
well  as  sample  inodes  of  a  particular  fd.  This  canbinaticn  allows  com¬ 
plete  mode  traces  for  all  values  of  sin©  and  fd  for  two  and  three  sub¬ 
stance  systems,  yielding  precise  and  reproducible  results. 

Due  to  repeated  use  of  the  subroutine  which  calculates  functional 
values  (on  the  order  of  100  to  500  times  per  data  point) ,  the  program 
should  be  run  on  a  computer  with  direct  access  to  the  central  proces¬ 
sing  unit  (such  as  the  DEC  PDP-MINC  11)  or  one  which  has  optional  use 
of  virtual  memory,  thus  cutting  confutation  costs  by  as  much  as  90%. 

By  actual  runs,  the  water/brass/water  systsn  results  of  Pitts  [8] 
have  been  confirmed.  Sane  of  the  calculated  results  of  the  air/ice/ 
water  system  have  been  verified  qualitatively  only.  Preliminary  lab- 
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oratory  runs  suggest  the  calculated  modes  are  correct.  Due  to  prob- 
lans  with  experimental  set-up,  determination  of  actual  parameters,  and 
recording  of  data,  insufficient  verification  has  resulted  fron  the  lab¬ 
oratory.  Figure  17  shows  the  theoretical  trace  of  the  average  Sato 
parameter  values  for  air/ice/water  and  two  experimental  points.  Hie 
frequency  used  was  1.84  MHz,  plate  thickness  was  approximately  1  irm 
(melting  of  ice  prevented  precise  measurement) ,  giving  an  fd  of  approx¬ 
imately  1.8.  Due  to  constraints  in  transducer  housing  and  set-up,  only 
angles  between  15°  and  60°  could  be  viewed.  TVo  non-specular  reflec¬ 
tions  were  observed,  at  26°  and  37°,  giving  the  two  points  shewn  in 
fig.  17.  Given  the  uncertainty  in  the  actual  parameters  of  the  systan 
and  the  experimental  problems,  the  points  correspond  well  with  the 
calculated  values. 

Future  research  should  include  continued  confirmation  of  calcu¬ 
lated  solutions  to  the  air/ice/water  system.  Methods  for  precise  de¬ 
termination  of  actual  parameter  values  such  as  shear  and  longitudinal 
w eve  velocities  should  be  pursued.  Though  Negishi  [22]  has  recently 
reported  on  negative  group  velocities,  this  area  should  be  investiga¬ 
ted  to  confirm  or  deny  the  existence  of  this  phenomenon.  Pitts  has 
recently  done  research  on  systsns  of  water/alumi nun  in  which  negative 
slopes  occur  in  the  mode  traces.  The  mere  presence  of  a  negative 
slope  may  or  may  not  indicate  the  existence  of  unusual  propagation 
characteristics . 

A  functional  basis  should  be  found  to  explain  discontinuous  inodes. 
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Though  no  physical  reason  is  apparent  for  such  discontinuities,  these 
inodes  may  be  the  "missing"  intercepts  representing  the  longitudinal 
wove. 

Should  the  modes  in  figs.  14  and  15  be  confirmed  by  experimenta¬ 
tion,  they  will  be  valuable  in  practical  applications  such  as  sonar, 
non-destructive  testing,  and  seismic  studies.  The  simplex  will  then 
have  proven  to  be  a  fast,  efficient  method  of  determining  solutions  to 
the  air/ioe/water  system,  as  well  as  any  other  similar  system  consist¬ 
ing  of  a  gas-solid-liquid  layered  structure. 
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A.  Simplex  FORTRAN  IV  acrputer  program  in  mode-trace  form.  Air/ 
Ice/Water  system  reflection  coefficient  for  simplex  subroutine. 
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*****  SIMPLEX  OPTIMIZING  PROGRAM  —  L I  NEAP  PROGRAMMING  « 

ORIGINAL  hGRK  CCNE  EY  FU/NG  AN C  ERCVER 
REVISED  BY  NC  ANC  F  1M6EPGLR  >  A t G  78 

FUHThER  KEVISEC  8Y  FIMEEFCER.  SEPT.  78-APkIL  79 

This  PCRM  CF  SIMPLEX  TRACES  A  MCCE  FROM  FC=G  TC  FCM.G 
OTHER  FORMS  SCAN  A  PARTICULAR  FC  FROM  SIN  7HE  |A*C  TC  1.0 

THE  VARIABLE  ARRAY  LFLAG  IS  SET  UP  TO  INCICATE  TFE  CRCER 
IN  mHICH  TFE  SIMPLEX  METHCO  REACHES  ITS  FINAL  ANSkER  — 

1.  A  REFLECT I CN  HAS  OCCLRREC 

2,  A  CCN TRACT ICN  FAS  CCCURREO 

3 1  A  NEW  MINIMUM  MAS  NOT  REACFEC 
4.  A  PARAMETER  REPLACEMENT 
S?  AN  EXPANSICN  GCCUKREC 
6.  A  FAILEC  CCNTRACT ION 

ThE  PC1NTS  AT  MHICH  THESE  ARE  PRINTED  IN  THE  PROGRAM  IS  A  GOOD 
INC1CATICN  OF  RhAT  TFE  VARIOUS  BRANCHES  ARE  CCING 

••••••«•«•*<««•••  A  SASIC  VERSION  CF  This  PROGRAM  HAS  ALSO 

titMtitiiiiiitii  SEEN  RR I  TEN  —  RANG  SPECIAL  FUNCTIONS  USEC 


Hill 


SEE  NELOER  ANC  MEAC  REFERENCE  IN  M.S.  THESIS 
OF  HIM8ERGER  FOR  DETAILS  CN  SIMPLEX 


DIMENSION  XI (2 J ,P ARM! 2. 3 1  .  PM  IN  43  )> C ENT  12  I 
DIMENSION  XRUItXC(2l>XE(2l 
DIMENSION  LFLAGISCOI 

««<  THE  ARRAY  SF  RETAINS  THE  LAST  ITTERATED  VALUE  TC 

«<«  CHECK  AGAINST  THE  PRESENT  CNE  —  FCK  TRACING  MODES 

DIMENSION  SPI5CJ 

INTEGER  K 

INTEGER  L 

CM 

F0*4.C5 

M-l 

SPINI«C. 19781 
F0«FC-C.C5 

IFIFC.LT. 1.751  GC  TO  5C0 

B-SPIkl 

GO  TO  3 

B«B*C.C1 

IFIB.GT.ISPI  IU0.1GI  I  GO  TO  1 
L-C 

*****  XK 11  IS  (HE  STARTING  REAL  SIN  GUESS 
*****  AKU  X1121  IS  THE  IMAGINARY  SIN  GUESS 
Xltl)«B 
X1I2  l«-0.01 

*****  ST  1  ANO  S 12  ARE  TFE  INCREMENTAL  VALUES  THE  SIMPLEX 

*****  USES  TC  FORM  IHE  • SIMPLEX*  ITSELF 

STl«O.COl 


nnnnnnn  r>  n  P>  n  r> 
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PARMtl. 2I-XU 1DST1 
PARMI2.1I-XK2I 
PARM(2,2I«X1(2) 

P  ARM (1.31= XI 11) 

PARM12»3I=X1C2I*ST£ 

5  COM  IME 

00  20  1*1.3 
OC  1C  J  =  1 ,2 
10  X1(J)*FA«M(J, 1) 

*****  THE  SUOKUUT  INE  IS  CALLED  UPCN  REFEATEOLY  1C  FIND 
*****  I HE  PUNCI 1CNAL  VALUE  CP  THE  VAR  1C IS  POINTS  IN  1HE  SIMPLEX 
CALL  FUNCT (XI. PI ,FD,N,£166) 

15  FMIN ( I )*F 1 
20  CONTINUE 
25  11*1 
12*1 
13*1 

OC  30  1*2.3 

IF,(FMIN(I1).LI.FHIN(1)I  11*1 
30  CONTINUE 
00  35  1*1.2 
35  CENT ( 1 1*0. 

OC  45  1*1.3 
IFd.EC.llI  GO  TC  45 
OC  4C  J* 1.2 

40  CENT(JI*CENT(J)*PAHM(J.I) 

IF  1FPIN(121.0T.FNIN(II»  CC  TC  41 
13*12 
12*1 

GC  TC  45 

41  13*1 
45  CONTINUE 

OG  50  1*1,2 

50  CENT!  I  l*CENT(I 1/2.0 

**********  TPE  •CRITER*  IS  THE  CRITERION  USEC  TC  CET ERMINE 
**********  THE  POINT  AT  hHICF  THE  SIMPLEX  HAS  I1ERAIEC 
****••*•*•  CLOSELY  ENOUGH  10  RHAT  HAS  BEEN  CALLED  ‘ZERO' 
CR1TER*1.0E— 06 
PCLE»C. 

UUIiMSi  THE  S  IMPLEX  ACTUALLY  ‘FALLS  UPON*  ITSELF 
411**11***  OR  COLLAPSES  —  HHEN  THE  DISTANCE  EETrEEN 
**********  THE  POINTS  OP  THE  'SIMPLEX*  IS  <  THE  CR1TER, 
*14*1*4***  YCU  HAVE  REACHED  A  ••ZERO".  NOTE  THAT  THIS 
*(********  *  •ZERO'  •  MAY  NOT  e£  IDENTICALLY  ECUAL  TC  0 
**********  IN  FACT,  IF  II  IS  APPROX.  ZERO.  IT  HAY  8E  CNLY 
**********  A  LCCAL  MINIMUM  —  NCT  PART  CF  THE  MODE  STRUCTURE 
00  60  1=1,3 

PGLE«PCLE+(FMIN(l I-FMIM  13 » 1**2 
60  CONTINUE 

POLE* SLR T ( POLE/ 2. C) 

I F (POLE. LE. CHI  TER  I  GC  TC  16C 
OC  £5  1=1,2 

XR(I  )  =  2.*CENT(  1  (-FARM  I,  ID 
85  Xim  =  XK(il 
L*L*1 

C  ******************  REFLECT ICN  ******************** 
LFLAG(L)=l 

CALL  FUNCUXl.FR  .FC.K.C166I 

IF(FM1N(  12). GE.FH.ANO.FR. GE.FMIM  13)1  OC  TC  130 
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1F(FR.lI.FNIN(13)I  GC  TO  140 

IFIFMlftl  id. GT.fH.ANC.FR. G1.FNINU2IJ  GC  TC  95 
OC  SC  1*1,2 

XCll  J*C.5APARM  1,  UI*0.5*CEN1 11) 

so  mu*»c(ii 

L*LM 

imuimunuit  ccntk  ac  r  icn  tiiuitniutuuit 

LFLAG <L)*2 

CALI.  FONCIIXl, FC. PC, K, 41661 
GO  TO  1J5 
95  L*L*1 

lflag!u*3***,,,4<  nc  mmtuuum 

00  ICO  1*1,2 

parm  i,u»=xmi) 

XC(  ll«.5»P4RM(  t,tll4.5«CEM(  I  > 

ioo  xim>icm 

CALL  FLiNCT(Xl«FC,FG,k,C166l 
105  IFIFR1NI  UI.GT.FC)  GC  TC  120 
L*L* 1 

titmtttmutiit  REPL.  PARAP.  tll*SSliJ*l|S***S« 

LFLAG  (0*4 
CO  115  1*1,3 
00  11C  J*l,2 

110  PARK  (J  ,1  )  *C.5«  I  PARM  ( J,  I )  *  FARM  |  -J  •  13  I  J 
115  CONTINUE 
GO  TC  5 

120  00  125  1*1,2 
X1(1I»CIII 
125  PAKN(I,UI=XC(l) 

FM1N1 1 1 )*FC 
GO  TC  155 
130  CC  135  1*1,2 
X1(II*XRUJ 
135  PARK  II.UIMM  IJ 
FKIN I  1 1 )  =  FR 
GO  TC  155 
140  CONTINUE 

OC  145  1*1,2 

XE  ( I  |*2,«XR  ( 1 l-CENT  ( 1 1 

145  xim*xem 

CALL  FUNCUX1,FE,FC,K,6166) 

IFIFPIM I3I.LE.FE)  GO  TO  156 
L-L4  1 

tlltlimilltllltA  EXPANSION 
Lf LAG (L 1*5 
00  150  1*1,2 

xim*xEm 

150  PARK!  I  ,  1 1 1 *XE ( 1 1 
PKINIIU-FE 

155  CONTINUE 
GC  TC  25 

156  L»L*1 


C 


C 

C 

C 


FAlLEC  CONTRACT ICN 


ittiiumu 


mumitmimt 
LPLAGlLI-6 
GO  TO  130 

l**M  THE  LFLAG  ARRAY  FAS  NOT  BEEN  FRINTEC  OUT  IN  THIS  FCRM 

*****  CF  THE  'IPFLEX  ROUTINE  IC  SAVE  l/C  TIKE 

««<  THIS  ROUTINE  CHECKS  TC  SEE  IF  INITIAL  2ERG  IS  AT  CORREC 
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C  ««<  START  INC-  PC  INI  —  AGAIN  fCO  MCCE  TRACING  »»> 

160  IF1C.NE.I.CI  GO  TG  161 

IFIPARMd,  m.LT.C. IflO.GR. PART  d,  I31.G1.C. 2201  GC  TO  16? 
h*2 

SP Ik l*P ARM ( 1 • X  3 1 

kRITEU.185)  (FARN(J,I3l,J*l,2l 

kRIT£(6,19C)  fMIM I3),*,f C.E 

C«Cd 

GC  TC  1 

C  <««  THIS  ACC  TINE  IS  USEC  FCR  CFECMNG  TC  SEE  IF  SANE  MQOE 

C  <<«<  HAS  BEEN  FCCNC  -  IF  NCT.  IT  OISCCLNIS  IT  ANC  TRIES  AGAIN 

161  CCNT1NCE 

IFiPARMIl,  I3J.GT.  ISP(V.)*0.Cl>.CR.PARMd,I3).LT.(Sf  Ikl-0.0311 
IGC  TG  162 

kR 1 TE (6 ,  18  51  (PARR (J,  131  , A*  1,2  I 
kR I T  E ( 6 , 19G 1  FMIM13),K,FC,B 

k*h»l 

SPIfcl-FARMd. 131 
C*C  Al 
GO  TO  1 

162  kA  ITE (6,200)  FC.F*«  IN d 3 1 , £ ,1  PARR  1 J,  13 1 , J*1 ,2) 

IF(G.EC.l.C)  GG  TC  2 

GG  TG  2 

166  1 F (G. EC. 1.0  I  GC  TC  169 

IFIPARMI 1, 13).GT.(SP(kl+C.Ql) .CR .P  ARM (1,I31.LT. I SF (k 1—0.03) 1 
I  GO  TO  16  2 

kRITE  (6,1931  FMIM  131 ,  FO ,  E  ,  I  PARR  ( J  ,  13 1 ,  J«l  ,2 1 
k-k*l 

SP(kl*PARH(l,13) 

C«G*1 

167  IF(C.EC.l.O)  GG  10  2 

168  GG  TG  1 

169  MR  I  TE (6  1 196  I 
GC  TG  160 

C  FORMAT  SECTION 

385  FORMAT  (///'  FINAt  FAAAPETER  VALLES  -  REAL  AND  IRAG.  PARTS  CF  SIN* 
I/I2E15.6II 

ISO  FORMAT  (/ ,  *  FUNCTION  VALLE**  ,j£l 5.6 ,  •  A  lTER*J.I7i*  FO’,F6.2, 

C*  B* *',£15.61 

193  FORMAT (/, •  NC  AFIER  3C0  I  TER .  •  .E1S.6 ,  •  *  FLNCT.  VALUE*, 

I/*  FD**  ,F6.2  ,  •  B*  '  , E 15.6 , /, 1  PARAM.  VALUES  KE ACFEC*  •  , / , ( 2E15.6 1 1 
196  FORMAT ( / ,  •  RINGING  FIRST  FOLE-NCT  CONVERGENT ,  CCNTINUING'I 

199  CONTINUE 

200  FORMAT!//*  PARAM.  VAl.  NCT  klTFIN  LIMITS  CF  PREVIOUS  VAl.*, 

I/,*  FC**,f6.2,*  FLNCT.  VAl.*  •  ,  E 15. t ,  •  E**,£15.6,/, 

I*  PARAMETER  VALLES*'  ,(2EL5.6I) 

310  CONTINUE 
500  STGP 
ENC 

SIBRCLT INE  FLNCT(X1,F1,FC,J<,*) 

C  SLBRCLT INE  FOR  NUMERATOR  CF  FLANG'S  REFLECT  ICN  COEFFICIENT 

C  FCR  MATER/1CE/A1R 

C  AAIAAAA  LGk  VALLES  FCR  VS  ANC  VL  LSEC  AAA AAV 

INTEGER  K 
DIMENSION  X 1 ( 2 1 

COMPLEX  A,C,EfG,R«C,R,S,T,Cl,ANG 
COMPLEX  STM, AST,  23,  21,  CC.Gl, GO. OCk 
COMPLEX  NUM 
COMPLEX  T1.T3 


U2*X1(2I 

AA,G«CPPIX(UI,02) 

0EMM.C2 

CEN3*O.OolX63 

DfcKP*C.S17 

VIM  •  5 

V3*C.34 

V 0*3.4 

VS*  1.5 

CXMC.COC,  t.CCCI 

P-3.1415S2t*FG/Vl 

AMVS/WXIMNG 

C*t 1. C-2.0* ( A**2 I »**i 

E*CSCRT  ((Vl/VCI*A2-(Af.  G**2)l 

G*CSCM(  (  Vl/VS»*»i-(  Afr£**2i  J 

HM.CCCM  (YS/Y1)**4>  *1  ANG**21*E*G 

G*CS1MP*EI 

r*ccc$xp*ei 

$*£SIMP*GJ 

T«CCGS(P*G) 

SVH*C*S*S»M*C»1 

ASr*c*c»r**-*K«5 

Y3*CSCRT<  (  V1/V3J**2-UXG**2I) 

YX*CSSRTU-IANG**i)> 

Z3*C1*(CEN3/CC,NPIME/Y3) 

ZX-ClMOENl/OENPIME/Vl) 

GC*Z1*K*T 

Gti*2  .CCO*  ISYH— 123*C*S)) 

GC«1C*S*T)+(MC*R) 

OQh*GLMAST*GCJ-(Zl-ZJJ*CC 

MM*GL*(ASr-GCI*(ZX-*Z3>*GC 

G?*REAMKtM) 

G8-A1NAGI.NUO 
FX*ABStG7)+AES<G8l 
K*K*  X 

1MK.XE.3DG)  Gt  U  1 0 
RETURN  X 
10  RETURN 
END 
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B.  Sanple  sinplex  output  for  node- trace  form.  Shown  are  points  that 
iterate  to  a  solution  to  the  reflection  coefficient,  points  that 
do  not  converge  within  a  set  ntanber  of  times  but  give  intermediate 
values,  and  points  that  do  not  belong  to  the  traced  mode  but  were 
given  during  the  course  of  the  iteration  for  informa¬ 

tion. 


L 


J 
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F I NAL  PARAMETER  VALUES  -  REAL 
(‘.318196E  CO  -C.  l41635E-0t 

FUNCTION  VALUE*  0.149012E-07  *  ITER*-  125  FD=  3.75 
8*  0.3264686  00 

NC  AFTER  3C0  ITEP.  0.13411CE-06  *  FUNCT.  VALUE 
FO*  3.70  B  =  0.318198E  00 

PARAM.  VALUES  REACHED* 

C.3L8662E  CO  -L  •  1285206-01 

NC  AFTER  300  ITER.  0.894070E-C7  *  FUNCT.  VALUE 
FO*  3.65  B*  0.308662E  OO 

PARAM.  VALUES  REACHED* 

0.297500.E  00  -Q.U3176E-01 


AND  IMAG.  PARTS  CF  SIN 


1 


FINAL  PARAMETER  VALUES  -  REAL  AND  IMAG.  PARTS  CF  SIN 
...  0. 284263 £  00  -0 . 959546E-C2 

FUNCTION  VALUE*  0.286847E-06  »  ITER*  96  FO*  3-60 

-B*  0.2975086  00 


.FINAL  PARAMETER  VALUES  -  REAL  AND  IMAG.  PARTS  CF  SIN 
C.268250E  tU  -0 . 7 74238E-C2 

FUNCTION  VALUE*  0.186265E-07  *  ITER*  109  FD*  3.55 
B*  0.284263E  00 

FINAL  PARAMETER  VALUES  -  REAL  AND  IMAG.  PARTS  CF  SIN 
0.240497E  00  -0. 563642E-02 

'FUNCTION  VALUE*  ’  0.447035E-07  #  ITER*  96  FD*  3.50 
S*  0.268258E  00 

J<C  AFTER  300  ITER.  0.175554E-06  *  FUNCT.  VALUE  _ 

FO*  3.45  B*  0.248497E  00 
PARAM.  VALUES  REACHED* 

_ 0. 223357E  CO  -0. 39866 IE-02  .  .  .  .  .  _ 


PARAM.  VAL.  NOT  bITHIN  LIMITS  OF  PREVIOUS  VAL. 

FD*  3.40  FUNCT.  VAL.*  0.S58754E-Q7  8*  0.223357E  00 

l PARAMETER  VALUES*  0.I89763E  00  -0.231284E-02  _l 
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C.  Sample  sinplex  output  for  scan  form  with  indicator  'flags'.  De¬ 
scriptions  of  flags  are  included  in  FORTRAN  IV  program,  appendix 
A. 


L 


_1 


FINAL  PARAMETER  VALUES  -  REM.  AND  IMAG*  PARTS  OF  SIN 
0.435409E  00  -0.I69694E  00 


NNOOOO 
—  **0000 
—  CM  O  Q  O  O 
m-(3000 
"■  I<1  O  O  O  o 
ifl-OOOO 
O  —  CM  O  O  O  O 
O  Ift  —  O  O  O  O 
X  M3  O  O  O 
m  ift  «  n  o  o  o 
O'  —  CM  o  o  o  o 
O'  m  —  o  000 
ffixNOOOO 
ff1  WhOOOO 
in  —  cm  o  o  o  o 
c>*  cft  —  o  ©  o  o 

•  "-00000 

o  in-o  o  00 

—  <M  a  o  O  O 
IflHOOOO 

—  cft  o  o  o  o 

n  ift  —  o  o  o  a 

CD  —  CM  O  O  o  o 

m  —  o  o  o  o 

OXNOOOO 

WM-OOOO 

•  —  CM  O  O  O  O 

-N-OOOO 

-4  CM  O  O  O  O 
CM  —  O  Cj  O  O 
1C  «*  CM  O  O  O  O 

a  'O  —  o  o  o  o 

U.  —  CM  O  O  Q  O 
CM  —  0  O  O  O 
CM  —  «0  O  O  O  O 
0«"0000 

—  —  cm  a  o  o  o 
*M  —  O  ©  O  O 
—  CM  O  O  O  O 
CM  -*  O  O  O  O 
—  CM  O  O  O  O 

fl  •««  —  <3  O  O  O 

oc  cm  cft  o  o  o  o 
Uf  —  -4  o  c.  o  o 
hMNOOOO 

—  —  -4  O  O  O  O 

«M  CM  O  O  O  O 

*«40UOO 
"4  CM  O  O  O  O 
l*»  *M  "4  o  O  O  O 
0-400000 
I  CM  *4  O  O  O  O 
Ul  *«4  CM  O  O  o  O 
—  CM  —  OOOOO 

—  —  CMOOOOO 

m  —  —  00000 

IMMNOOOOO 

00"4"4000o0 

—  —  OOOOO 

•  mcMooooo 

O  —  —  OOOOO 

eft  cft  o  o  o  o  o 

. 4-NOUOO 

i«  cm  —  ©  o  o  o 

0  —  —  ♦  o  O  O  O 
UJIflOlMOOOO 
3  —  —  —  oooo 
jino-oooo 

<•4*4^0000 

>cocncMoooo 

•4  *4  —  o  o  o  o 

zn-ooooo 

O-iAmOOOO 
mCC4oOOOO 
Jb—  —  —  OQOO 
IJlAlMNOOOO 
Z———  OOO© 
3«CMCMOOOO 
,11444  OOOO 


UJ 

z 

— 

eft 

u. 

3 

0 

-J 

•« 

(ft 

> 

H- 

• 

a 

< 

b~ 

b. 

U 

z 

• 

3 

0 

UL 

4 

n 

< 

£ 

** 

>6 

0 

0 

z 

1 

< 

UJ 

m 

0 

_J  0 

<r 

0 

*t  0 

CM  O 

UJ 

f4  O 

UJ 

a  uj 

CM 

* 

<♦ 

CM  CU 

O' 

1  O' 

*  O' 

•o 

vO 

O'  Ci  >0 
C*  Uj  — 
*«  X  • 
•MJO 

k  *n  I 

uj  o  uj 

a. 

•4  O 

c/>  o 
O  H  UJ 
O  03  3  UJ 

O  JO 
—  O  <  — 
eft  >  •* 

•  •  Ift 

u.  z  + 

<  <t  • 
n  oc  O 
o  o  < 
z  u.  a 


CM  O  O  o  o  o 
•400000 
<N  O  O  O  O  o 
<400000 
CM  O  o  o  o  O 
•400000 
O  <M  O  O  O  O  o 

o  <4  o  o  o  0  o 

cm  o  o  e  o  o 

UJ  -4  O  C  •  O  C3  o 

O'  CM  O  O  O  O  O 
<7*  *4  0  Q  O  O  o 
O'  CM  O  O  O  O  o 
O'  <4  O  o  o  o  o 

N  O  O  O  O  O 
H  *4  o  O  000 

•VOODOO 

o  “•  o  o  o  o  0 

-400000 
<0  —  0000 

—  >r  o  o  o  a 
U  «NOOOO 
0-4-40000 

♦  o  o  o  o 
o  *4  CM  O  o  o  o 
m  4  -  o  o  o  o 

•  in  cm  o  o  o  o 

———OOOO 
Cft  cm  o  o  o  o 

—  —  o  o  o  o 
n  »m  ♦  o  o  o  o 
a  —  cm  o  o  o  o 
U.  <0  -4  O  O  O  o 

■400000 
O  CM  -4  O  O  O  O 

Cft  -4  ♦OOOO 

•4  CM  CM  O  O  O  o 
•4  —  0000 
tft  CM  O  O  O  O 
-4  —  O  C  O  o 
CM  CM  o  O  O  O 
II  —  —  0000 
KflOOOOO 
tu  —  —  o  o  o  o 

—  cm  o  o  o  a 

—  •0  —  0000 

—  CM  O  OOO 

*0—0000 

—  CM  O  O  O  o 
«0  cm  —  O  o  o  o 
O  —  CM  o  O  O  o 
|  CM  —  O  O  O  o 
UJ  —  <0  o  o  o  o 
•0  cm  —  OOOOO 

—  —  INIOOOOO 
43  cm  —  OOOOO 
»  —  CMOOOOO 
O  —  —  OOOOO 
CMIftcftOOOOO 

•  —  —OOOOO 
OlftCMOOOOO 

•4—00000 
Cft  cm  O  O  O  O  O 

——OOOOO 

n  eft  <m  o  o  o  o  o 

UJ  —  —  OOOOO 
3lft  —  OOOOO 
«J  —  CMOOOOO 
<  Cft  —  OOOOO 
>  —  CMOOOOO 
Cft  —  O  O  o  o  o 
Z  —  CMOOOOO 
O  —  —  OOOOO 

—  CftCMOOOOO 
K  —  —  OOOOO 

ucnmooooo 

Z  —  —  OOOOO 
OlftCMOOOOO 
IL  —  —  OOOOO 
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